% fig3a.m ——— 绘制 Fig.3(a1–a3)
clear; clc; close all;

%% —— 模型参数 —— 
a     = 0.5;
b     = 0.215;
beta  = 0.2;
A     = 0.95;
omega = 0.61;
alpha = 1.7;

%% —— 仿真设置 —— 
dt      = 0.01;             
T_total = 1000;      % 总时间
T0      = 200;       % 丢弃瞬态时间
tspan   = 0:dt:T_total;
X0      = [0.01; 0.02; 0.02];

%% —— 积分 —— 
f = @(t,X) rhs(t, X, a, b, beta, A, omega, alpha);
X = rk4(f, tspan, X0);

%% —— 丢弃前瞬态 —— 
idx = tspan >= T0;
t_plot = tspan(idx);
x_plot = X(1, idx);
z_plot = X(3, idx);

% —— 计算能量 —— 
H_plot = energy(X(:,idx), alpha, a, b);
H_avg  = mean(H_plot);

%% —— 绘图 —— 
figure;

% (a1) x–z 相图
subplot(3,1,1);
plot(x_plot, z_plot, 'k', 'LineWidth',1);
xlabel('x'); ylabel('z');
title('(a1)','FontSize',12);
axis([-1.2 1.2 -2.5 2.5]);
box on;

% (a2) x 随时间 τ 演化
subplot(3,1,2);
plot(t_plot, x_plot, 'r');
xlabel('\tau'); ylabel('x');
title('(a2)','FontSize',12);
xlim([T0 T_total]);

% (a3) H 随时间 τ 演化
subplot(3,1,3);
plot(t_plot, H_plot, 'b');
xlabel('\tau'); ylabel('H');
title(sprintf('(a3) \\langleH\\rangle=%.4g', H_avg),'FontSize',12);
xlim([T0 T_total]);


